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Abstract. In this lecture I discuss the bulk surface heterogeneity of rotating stars, 
namely gravity darkening. I especially detail the derivation of the ca-model of Espinosa 
Lara & Rieutord (201 1|) , which gives the gravity darkening in early-type stars. I also 


discuss the problem of deriving gravity darkening in stars owning a convective envelope 
and in those that are members of a binary system. 


1 Introduction 


As it has been much discussed in this school, surface inhomogeneities of stars 
are more and more frequently detected due to the increasing sensistivity of the 
instruments. If correctly understood, and therefore modeled, these data may 
open new windows on the interior or on the history of stars. 

The purpose of this lecture is to first briefly review the processes that lead 
to such inhomogeneities and then to focus on a very fundamental one, namely 
rotation. 

In the ancient times, when the eye was the only optical instrument to observe 
Nature, the sun was thought as a pure uniform bright disc. The invention of the 
telescope by Galileo ruined this idea, showing that the sun was spherical with 
spots on it. Presently, it suffices to have a look at images of the sun taken at 
short wavelengths to understand that its surface is certainly not uniform. Such 
images actually reveal that the magnetic fields are a prominent cause of this non¬ 
uniformity. It looks like a mess which even impacts the distribution of surface 
flux and temperature. A closer look at the magnetic structures but also below, at 
the photospheric level, shows that all these heterogeneities evolve with time. On 
the photosphere, turbulent convection features the surface with two important 


scales: granulation and supergranulation (Rieutord & Rincon 2010). Even the 
bulk surface rotation is not uniform. This differential rotation, known since the 
nineteenth century, with fast equatorial regions and slow polar regions, is now 
understood as driven by Reynolds stresses coming from the turbulence in the 
solar convection zone. 

Hence, the surface of the sun teaches us that we should expect non-uniform 
velocities, temperatures, flux and magnetic fields at the surface of all low-mass 
stars. But one consequence of the strong mixing imposed by turbulent convection 
and the ever changing magnetic fields is that the solar photosphere has a uniform 
chemical composition! 
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But uniform chemical composition is certainly not possible when turbulent 
convection disappears and no longer mixes the surface layers, that is when we 
consider stars of higher mass with an outer radiative envelope. There, com¬ 
bination of magnetic fields with microscopic diffusion processes (gravitational 
settling or radiative acceleration) may on the contrary lead to chemical surface 
inhomogeneities ( Vauclair fe Vauclair|1982 Alecian|2013 Korhonen et al.|2013 ). 
But even when magnetic fields are absent at their surface, early-type stars are 
still endowed with a non-uniform surface: absence of magnetic field is indeed 
correlated with fast rotation, a feature that makes the polar caps brighter than 
the equatorial regions. 

From the foregoing presentation we see that three processes, intrinsic to each 
star may lead to surface inhomogeneities: rotation, convection and magnetic 
fields. There is a fourth one, but extrinsic to the star itself, namely binarity. 
A companion indeed raises tides, illuminates one side of the star or may even 
transfer mass. 

Within these four physical processes that make the surface of stars not uni¬ 
form, we shall concentrate on the most simple, namely rotation, which, a priori 
leads surface variations that only depend on the latitude. We shall discuss in 


detail this very basic physical effect, leaning on the recent works of Espinosa 


Lara & Rieutord (2011 2012). 


2 The energy flux in radiative envelopes of rotating stars 

2.1 von Zeipel 1924 


In a seminal paper, von Zeipel (1924) showed that a rotating star may be brighter 


at the poles than at the equator. This result is quite simple to derive if we 
assume that the star is barotropic and that the energy flux is given by Fourier’s 
law. Indeed, if the star is barotropic, meaning that its equation of state can be 
simplied to 


P = P(p), 

it implies that there exist a hydrostatic solution in the rotating frame and that all 
thermodynamic quantities can be expressed as a function of the total potential 
(i.e. gravitational plus centrifugal). Hence, one writes 

P = p($), T=T($), etc. 

Then, using Fourier’s law to derive the heat flux, one finds 

Frad = -xVT = -xo^ro?) = K($)g cS 
whence von Zeipel law 


TeS 


Kg'Jt 


on the surface = Cst 
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This result is simple but incorrect. Indeed, barotropic stars are realized in two 
cases: either the star is isentropic and thus fully convective or it is isothermal but 
this can hardly be the caseQ So the closest case that may match the barotropic 
state is that of a fully convective star, but in such a case the flux cannot be 
derived from the Fourier’s law. In fact, these hypothesis (barotropicity and heat 
diffusion) lead to a contradiction. Indeed, we may note that the total potential 
and the temperature would verify in the envelope of the star, 

f Div(yVT) = 0 

( A<P = 4nGp + 21? 2 y > 

which leads to 


Di v(x(<2>)T'(<I>)V<Z>) = 0 


AirGp + 21? 2 + (\n(xT'))'g 2 eS = 0 


On an equipotential, p is constant as well as (ln(xT / )) / , but g e g is not constant. 
Hence, the latter equation is impossible. The reason for that is that for a rotating 
star where heat is transported by diffusion, a barotropic state cannot be and 
should be replaced by a baroclinic state. In such a state, isobars, isotherms 
or equipotential are all different, not very different, but different. This is the 
normal state that comes from the fact (basically) that temperature, pressure, 
and gravitational potential all obeys different and independent equations. The 
barotropic state is therefore rather peculiar (but see Rieutord [2006 , for a more 
detailed presentation). 

Now one may wonder if it is possible to derive the dependence of the flux 
with latitude for a rotating star without computing the whole stellar structure 


and the associated flows as in Espinosa Lara & Rieutord (2013). Fortunately, 


this is indeed possible as Espinosa Lara & Rieutord (2011) have shown. It is 


not as simple as von Zeipel law, but it has the merit of relying on controllable 
hypothesis. 


2.2 The idea of the w-model 

In the following we shall first restrict ourselves to the case of early-type stars, 
that is to stars that have a radiative envelope around a convective core. We’ll 
discuss the case of convective envelope in the next section. 

Within the envelope of a star the flux just obeys: 


DivF = 0 (2) 

namely energy is conserved and there are no energy sources. 

This is a single equation, not enough to determine the two components, 
(. F r ,Fg ), of the flux, but if we add a constraint to the flux we may find it. We 
thus assume that the flux is anti-parallel to the effective gravity 

1 Some models of prestellar core use this hypothesis, sometimes. 
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Fig. 1. Misalignment between pressure gradient and flux for a configuration with a 
flatness ~ 30% (from Espinosa Lara fe Rieutord|2011 1. 


F = -f(r,9)g cB (3) 

In order to avoid an additional unknown, we shall take the effective gravity g eB 
as given by the Roche model. In such a case we shall see that the flux function 
/ can be determined and that the latitude variation of the flux depends on a 
single parameter u> defined as the ratio of the angular velocity to the keplerian 
angular velocity at the equator. In other words, the flux depends on 



Thus, we shall call this model the u) — model to emphasize the crucial role of the 
reduced angular velocity to. 

However, before going any further, we may wonder whether the assumptions 
are strong or not, especially ([3|. 

In a radiative zone, the configuration is baroclinic so vectors are surely not 
aligned but fortunately we can now revert to 2D-models to get an idea of the 
misalignment. As shown in Fig. |TJ the misalignement remains small, less than a 
degree, even if the star rotates close to criticality. 

Thus, even if the envelope is the seat of baroclinic flows, the misalignment 
is small. Actually, the baroclinic torque (VP x Vp)/p 2 does not need a strong 
misalignment of the vectors to be efficient at driving baroclinic flows because 
the two gradients (of pressure and density) are already quite strong. 

Let us pursue somewhat. From and ([3]), we have 
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Fig. 2. Left: Density profile of a M=5 Mq ZAMS, non rotating star (solid line) together 
with that of a n=3/2 and n=3 polytropes. Right: The interior gravity of the same 
stellar model (solid line) together with the — 1/r 2 Roche model (dashed line). The 
ZAMS model is an ESTER model (e.g. Espinosa Lara fc Rieutord|2013|l 


DivF = 0 Div(/W) = 0 

thus 


Soft • V In / = -2I2 2 

because A<f> = — 2L? 2 ; hence, 


(5) 


a In/ _ 21? 2 

d£ g e s 


( 6 ) 


where we introduced the local vertical coordinate £. Equality © shows that 
, and therefore /, has latitudinal variations similar to those of g e g. Hence 
the horizontal variations of the flux cannot be given by von Zeipel law. In other 
words T e ff/glff 4 cannot be constant. 

The second hypothesis is the use of the Roche model. This model assumes 
that the whole mass of the star is concentrated at the centre thus leading to 
a gravitational potential in 1/r 2 everywhere. For the regions we are interested 
in, namely the envelope of early-type stars, this is a rather good approximation 
since these stars are usually said to be “centrally condensed”. In Fig. [2] we show 
the density profile of a non-rotating 5 Mq star along with the profile of two 
polytropes. We see that the n=3 polytrope represents fairly well the density 
profile of the star and that the n=3/2-polytrope, which is a very good model for 
fully convective stars, is much less “centrally condensed”. Hence, gravity in the 
outer envelope of an early-type star is well represented by the Roche model (see 
Fig. fright). The interior discrepancy with the Roche model has no consequence 
for our purpose. 
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2.3 The derivation of f(r,0) 

f is given by ([5| but we first need to scale this function so as to introduce a 
non-dimensional function F that accounts for the radial and latitudinal variation 
of the flux. This is easily done if we observe that near the star’s centre 


So that we may set 





GM 


with 


firi 9) 


L 

A-kGM 


F(r, 9) 


(7) 


lim F(r, 9) = 1 


r —>-0 


Then, we scale the gravity with GM/R/f and the length scale with the equatorial 
radius R e . The scaled angular velocity is therefore given by 



At this point we should underline that the angular velocity is scaled by the 
keplerian angular velocity given by the equatorial radius. It is often the case in 
the literature that the scale of angular velocity is the critical velocity associated 


with the Roche model of the considered mass M (e.g. Monnier et al. 2012 for 


instance). This gives a different ui (i.e. fraction of critical velocity). We give in 
appendix the relation between these two ways of appreciating angular velocity. 

We now proceed to the derivation of F the scaled version of /. From 0 and 
<[5j) we get 


— r sin 


dF dF 

— -sin(9cos6> — = 2 F 
or 09 


With F(0,9) = 1 we have all the elements for solving the equation for F. 
First, we solve for lnF 1 , namely, 


1 


— r sin 2 9 


9 In F . „ d In F „ 

,,-sin 9 cos 9 ——— = 2 . 

or 09 


If we set In F = In G+A(0), so that A(9) removes the RHS of ( 
find that 


( 8 ) 

we immediately 


A'(9) = —2/ sin 9 cos 9 => A{0) = — ln(tan 2 9) . 
But we still have to solve the homogeneous equation, namely 
1 . , 


— r sin 


51nG . 5 In G 

— sin 9 cos 9 ——— = 0 


dr 


09 


(9) 
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Such a first order partial differential equations is solved by the method of char¬ 
acteristics. We therefore look for places where InG is constant. These places are 
called the characteristics curves of G. They are such that 


9InG , 9InG „ 

dr T+ 86 d 

but G also verifies M so that we can eliminate dl g G and dl g g G and get 


- — r sin 2 8J dO + sm8cos ddr = 0 (10) 

which is the equation of characteristics. 

We first observe that we may multiply this equation by any function H(r, 0) 
without changing anything. So we may also look for h such that 


dh 

dr 


H sin 8 cos 8 


dd 


— r sm 


(ii) 


where H needs to be chosen so that this system can be integrated. After trial 
and error, we find that H = ui 2 r 2 cos 8 cot 8 is the right function. Thus 


and the solution is 


dh 

dr 


= ui 2 r 2 cos 3 8 


dh 

88 


cos 2 8 
sin# 


w 2 r 3 cos 2 9 sin 8 


( 12 ) 


h(r,8) = ^w 2 r 3 cos 3 8 + cos 6 + lntan(0/2) 

The curves h{r, 8) = Cst are the characteristics. Note that the polar equation of 
a characteristic, namely the dependence r = r(9), is just implicitly known, and 
depends on the chosen constant. 

Now, we know that In G or G is constant on the curves where h(r, 8) = Cst. 
So we can write 


G = G[h) . (13) 

It means that the variations of G with (r, 8) are through those of h(r, 8) only. 
So we find that 

InF = lnG(/i) — In tan 2 8 or F = ^’ —— 

v ' tan 2 8 

This is the solution of the partial differential equation, but it is up to an arbitrary 
function G(h) that we should determine. For that, we need to revert to the 
boundary conditions, namely that F(0, 8) = 1. We thus need to impose 










Michel Rieutord 


or 


G(h(0,0)) 
tan 2 8 


(14) 


G(cos0 + lntan(#/2)) = tan 2 0 (15) 

for all 9. This is certainly a weird expression of G, but actually it is sufficient. 
Let’s introduce the function h 0 such that 


Hence, we have 


or 


h 0 {8) = cos 9 + lntan(#/2) 


(G o h 0 )(9 ) = tan 2 9 


(16) 


G o ho = tan" 
so formally, the solution for G is 


G = tan 2 oh Q 1 


G(r, 9) = tan 2 (h 0 1 (h(r, 0))) 
To make it more understandable, we set 


so that 


or 


it> = h 0 ( h(r,9 )) 


ho (ip) = -w 2 r 3 cos 3 9 + cos 9 + lntan(0/2) 

O 


(17) 


cosip + lntan(i/’/2) = -w 2 r 3 cos 3 9 + cos9 + lntan(0/2) 

O 


(18) 


which is a transcendental equation for ip. However, it is not difficult to solve 
numerically (we know that when r or w are small ip ~ 9). So finally we find 


F(r, 0) = 


tan 2 {ip{r, 9)) 
tan 2 9 


(19) 


where ip(r, 9) is given by 


2.4 Two interesting latitudes 

F seems to be singular at the pole (9 = 0) and at the equator (9 = 7r/2). Let us 
explore these two latitudes. 


Starting with the pole, we see that if 9 <C 1, then, from (18), we find that 
ip 1 as well. Indeed, for small values of the angles we have 
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1 + lntan(-0/2) ~ r w 2 r 3 + 1 + lntan(#/2) 

3 

so that 

ip ~ 6»e w2r3/3 (20) 

and 

F(r, 0) = e 2 " 2r3/3 (21) 

which gives the values of F along the rotation axis. 

The equator is more complicated. We need to know that if £ <C 1 then 

In ^tan ^ — ej = — e — ^e 3 — • • • 

With this asymptotic expansion we find 

F(r, tt/ 2) = (1 - w 2 r 3 )" 2/3 


2.5 The final solution of the o>model 

Back to the definitions we started with, we can express the flux with the effective 
gravity in the following way: 


F = -^GM Fi “’ r ’ e)g °« 

so that we also get the effective temperature 


( 22 ) 


T e fi 


( L \ 1/4 / tan ip 1/4 

V4t tcjGMJ V tan 9 9eS 


(23) 


From this expression, we see that the function \/tan ip/ tan 9 shows the deviation 
from the von Zeipel law. 

Noting that 


GM ( e r 2 • n \ 

Sm 0e s ) 

for the Roche model (e a is the unit radial vector of cylindrical coordinates and 
e r that of spherical coordinates). We find that the flux is given by 


F = -^(-^+^rmn9e.)F( Ut r t 9) (24) 

which shows that it depends only on u> and a scaling factor , hence the 
name “w-model”. 
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Fig. 3. Left: Scaled effective temperature as a function of colatitude for a 3 M@ model 
at 17 = 0.9 Qk- The solid line shows the prediction of the simplified model, ’pluses’ show 


the prediction of a fully 2D ESTER model including differential rotation (see Espinosa 
Lara & Rieutord 2013| for details), while the dashed line is for the von Zeipel law. 


Right: With the same symbols as on left, the effective temperature as a function of the 
effective gravity (plots from Espinosa Lara fe Rieutord|2011 1. 


2.6 Comparison with 2D models: a test of the (3- and uj- models 


After the foregoing mathematical developments we certainly would like to com¬ 
pare the results of this modeling to more elaborated models. For this purpose, 
we compared the latitude variations of the flux with the prediction of fully two¬ 


dimensional ESTER models (Espinosa Lara & Rieutord 2013). We recall that 


ESTER models give a full solution of the internal structure of a rotating early- 
type star including the differential rotation and the meridional circulation driven 
by the baroclinicity of the envelope. They also include the full nricropliysics 
(opacity and equation of state) from OPAL tables. Fig. [3] and [4] show that the 
w-model matches very well the output of the full ESTER models. Moreover, we 
also note that the dependence of the effective temperature versus gravity is close 
to but not exactly a power law. 

Observational data often show the polar-equator contrast in effective tem¬ 
perature in terms of the exponent /3 defined as 


T eS (X <4 (25) 

We shall call this approximate modeling the “/3-model”. Actually, note that ( f25j ) 
demands that 


P 


d In T eff 


d In g eS 


r=R(8) 


(26) 


where R(9) is the radius of the star at colatitude 9. Since the relation between 
T e gf and <7 e ff is not a power law, /? is not constant on the surface of a rotating 
star. It varies between two extreme values that we can also compute. 
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Fig. 4. Variation of the ratio of effective temperature at pole and equator as a function 
of the flatness of the star. Symbols are the same as in Fig. [ 3 ] (plots from Espinosa Lara 
|fe Rieutord||20lT| ). 


To make things simpler we therefore define the 6-exponent as follows: 


T = T 



ln(T e /T p ) 

In (g e /g P ) 


(27) 


where the indices e and p refer to the equator and pole respectively. T and g 
designate the effective temperature and effective surface gravity. 

From the polar and equatorial expression of the flux, we get 


F e = (1 - w 2 ) 2/3 g e and F p = e 2u,2r p /3 g p 
for the w-model, while, from the Roche model, 

— =r 2 (l-w 2 ) with r p = 1 
g p 1 + w 2 /2 

where r p is the polar radius. So we find 

(TeV_ (1 ~~ ta 2 ) 1 / 3 -2u 2 r!/3 

\T p ) (1 + W 2 /2) 2 


,1 1 ln(l — w 2 ) + uj 2 r p 

4 6 ln(l — w 2 ) — 2 ln(l + w 2 /2) 1 ; 

We plotted in Fig. [5] the values of 6 with increasing values of the flatness 
(namely with increasing rotation). In this figure we see that the 6-exponent is 
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0.00 0.05 0.10 0.15 0.20 0.25 0.30 

epsilon 


Fig. 5. The /3-values from various models: The solid line shows the 6-exponent of the 
w-model, while the dot-dashed line shows the corresponding ESTER model. The extra 
dashed lines give the range of /3 values spawned at the stellar surface by a /3-model. Data 


from interferometric observations of some early-type stars are shown (from Domiciano 
de Souza et al.|2014 l. 


close to a linear dependence b = \ — \e up to £ = 0.3. However, note that since 


the true dependence is not a power law, /3, as given by (26), varies at the surface 


of a given star. We also show in Fig. [5] its range of variation and it is clearly 
not negligible when e is larger than ~ 0.15. This means that if we had access to 
a very high spatial resolution of the stellar surface we would find different /3’s 
whether we look at the pole (large values) or at the equator (low values). 

If the /3-model is a poor representation of the latitudinal variation of the flux, 
can we devise a better one? Surely, a decomposition of the effective temperature 
on the spherical harmonics basis, namely 

r e (0) = £4yr 

Z,ra 

has the advantage of being model independent. The coefficients of the expan¬ 
sion are the results of observations. Such an expansion is already used for the 
description of spotted stars for the reconstruction of their magnetic fields (e.g. 
Do nati et al.| [2006 but see also the lecture of O. Kochukhov in this volume). 

In Fig.[5j we also show the observationally derived values for a few early-type 
stars observed with interferometers. The matching is quite remarkable, even if 
some cases like a Cep certainly need a more detailed study. 

To finish with the case of early-type stars, let us consider the case of small 
rotations. We may first derive the linear dependence of the 6-exponent with e. 


From (28) we get 
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, 1 1 2 
b = - — -CO 
4 6 


e>(w 4 ) 


or 


&=---£■ 
4 3 


0 ( 8 2 ) 


(29) 


where we observed that e = 1 —r p . This expression shows that in the limit of small 
rotation we recover von Zeipel law. To understand the origin of this property, it 
is useful to reconsider the to -model and the expression of the function F(r,9). 
Let us first solve (18) in the limit uCl. This yields 


i/j = 6 + 1-uj 2 sindcos 0 + O(to 4 ) (30) 

O 

From this relation, we derive the asymptotic expression of F(r,9) at low w, 
namely 


F(r,9) = 1 + ^w 2 r 3 
o 

The latitudinal dependence has disappeared. Hence the latitudinal variations 
of the flux are those of the effective gravity. Therefore, von Zeipel law applies 
at low rotation rates. We can understand this result, if we recall that in the 
limit of zero rotation, the star is spherical and all surfaces of constant pressure, 
temperature, etc. are spheres so that we can consider the gravitational potential 
or the pressure as the independent variable. Thus we recover a kind of barotropic 
situation where one can use a relation between pressure and density, and derive 
a von Zeipel law. 


3 The case of convective envelopes 


3.1 Lucy’s problem 

In the sixties it was realized that gravity darkening was very important for the 
interpretation of light curves of contact binaries (like the W UMa-type stars). 
But most of these stars are low-mass stars, thus with a convective envelope. The 
use of von Zeipel law, which is based on heat diffusion, was therefore doubtful. 
So Lucy (1967) asked: “What is the gravity-darkening law appropriate for 


late-type stars whose subphotospheric layers are convective?” Lucy’s reasoning 
was the following. 

In the convective envelope of a rotating star, if we go deep enough, we should 
reach a medium of constant entropy. This value should be the same whatever 
the latitude. ID models show that the entropy jumps from a minimum near the 
surface (where the convective driving ceases) to a plateau in the deep layers 
where convective mixing is efficient (see Fig. [6]). Lucy argues that the value of 
the entropy s on this plateau is a function of the surface gravity g s and effective 
temperature T e ff. He thus writes 


s = s(g s ,T e ) 

In the case of a rotating star, where g s and T e g vary, we must have 


( 31 ) 
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STEIN & NORDLUND 
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Fig. 6. Thermodynamic profile of the Sun according to Stein & Nordlund 1998. 


s{g s ,T e ) = s 0 (32) 

where sq is the entropy on the plateau. If we differentiate this expression with 
respect to g s and T e , we find that 


ds_ 

dg s 


dg s + -§fr dT e = 0 


in the deep layers of the rotating star. Since we admit that T e oc </f, then we 
have 


ds ds 

d In g s + d In T e 


(33) 


Thus, if we are able to evaluate the values of the above partial derivatives, we 
can obtain f3. For that, Lucy considered various ID neighbouring models (we do 
not know how the variations were made), and evaluated the partial derivatives 
so as to find /3. Using five stellar models (3 with M=1M 0 , 2 with M=1.26M 0 ), 
he found that 


0.069 </3< 0.088 
Lucy adopted j3 = 0.08 as a representative value. 

3.2 A new derivation of Lucy’s result 

It is interesting to note that Lucy’s results may be derived from simple consid¬ 
erations on one dimensional stellar models in the solar mass range. 
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Let us first recall that the surface of a star is usually determined by a surface 
pressure given by 

P = 3 ^ ^ 

where g s is the surface gravity and k an average opacity. This boundary condition 
comes from the assumption of hydrostatic equilibrium of the atmosphere, namely 

dP __ 1 dP _ g dP _ g 

dz ^ '' pn dz k ^ ** dr k 

where the last relation is integrated from the zero optical depth down to r = 2/3. 
In the range of density and temperatures typical of the solar type stars, opacity 
may be approximated by a power law of the form: 


k = KqP^T 


(35) 


For instance Christensen-Dalsgaard uses p = 0.408 and s = —9.283 for the sun 
(e.g. Christensen-Dalsgaard fc Reiter| [l9951. 

Now in convective envelopes, the variation of pressure and density are related 
to temperature through 


P oc T n+1 and p oc T" . 

namely with a polytropic law with n = 3/2. 

Using the foregoing power laws for the opacity, pressure and density, we can 
express gravity as a function of temperature. We find that 

g oc T n ^ +1)+1 ~ s 

Identifying temperature and effective temperature, we find a gravity darkening 
exponent which reads: 


P 


1 

n(p + 1) + 1 — s 


(36) 


Using Christensen-Dalsgaard’s solar values and n = 3/2, the foregoing expression 
yields 


P ~ 0.0807 

which is precisely the value found by Lucy. This is no surprise since Lucy used 
models similar to solar models, so the power law fit of Christensen-Dalsgaard is 
appropriate. 

This derivation clearly shows that this /3-exponent, as defined by (33), de¬ 
pends on the chemical properties of the surface through the opacities. 
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3.3 Can Lucy’s law represent a gravity darkening effect? 


The foregoing derivation of Lucy’s result enlights us on the origin of Lucy’s value 
of the ft exponent. We see that it is essentially due to the strong dependence 
of opacity with temperature in the surface layers. Since the values for /r and s 
are chosen to fit the table values in some range of density and temperature, we 
understand that Lucy’s result applies only to stars similar to the Sun, in terms of 
gravity and effective temperature. We may note, as [Espinosa Lara & Rieutord 
(2012), that if the opacity law extend in the deep layers so as to control the 


structure of the envelope and leads to a radiative one, then /3 = 1/4 because 
the polytropic index is n = (s + 3)/(p, + 1). We recover the previous result for 
non rotating radiative envelopes. We see that when the opacity is such that 
the polytropic index is less than 3/2, and the envelope is convective, the /3 is 
governed by the opacity of the surface layers, those which are assumed to be 
transparent and fixing the atmosphere. The structure of the envelope is close to 
the adiabatic index n=3/2. 

Now we wish considering the case of rotating stars. The question is whether 
we can use the foregoing value of the exponent, if we consider a fast rotating 
star of solar type. A first obstacle is the validity of the boundary condition 


(34), which relies on a hydrostatic equilibrium. When rotation is present such 


an equilibrium is impossible because of baroclinicity (for the same reason as 
the origin of the so-called von Zeipel paradox, see [Rieutord 2006). The proper 


boundary condition, replacing (341 should be derived from 


v ■ 'Vv = —VP — V<£ 

P 

where v is the fluid velocity in an inertial frame. Basically the flow is a differential 
rotation plus some weak meridional currents. The important point is that the 
differential rotation is latitude dependent. Hence, if we were to use some pressure 


boundary condition like (34), we should expect some extra variations from this 


latitudinal differential rotation. But this is likely not the whole story as we shall 
discuss it now. 

If we consider the deep convective envelope of a rapidly rotating star, we 
might consider too contradicting effects. First the Coriolis effect: analysis of 
a linear stability of a convectively unstable layer shows that polar regions are 
less unstable than equatorial ones. This a consequence of the presence of the 
Coriolis force. This force indeed prevents variations of the velocity field along the 
rotation axis (the so-called Taylor-Proudman theorem). It shows up in numerical 
simulations as convective rolls parallel to the rotation axis near the equatorial 
regions (Glatzmaier & Olson 1993). For stars this may imply that the convective 
flux is larger in the equatorial plane than in the polar region, thus meaning a 
negative /3. However, in the equatorial plane the effective gravity is less, and so is 
the buoyancy force. This is the effect of centrifugal force, which therefore points 
to more flux in the polar region (thus for a positive ft). The conclusion of the 
foregoing argument is that nothing is clear. We may only guess that if Lucy’s 
law applies, it is for slowly rotating stars of solar type only. This is a deceptive 
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Fig. 7. Schematic representation of the primary star with filling factor p = 0.8. The 
position of the Lagrange point L\ is shown on the 2 -axis that joins the centre of the 
two stars (from Espinosa Lara fc Rieutord|2012 1. 


conclusion since we may have interesting data only on fast rotators or evolved 
stars with weak self-gravity. In addition, the previous remarks do not mention 
the magnetic fields that are almost unavoidable in late-type stars. 


4 Binary stars 


Binary stars is another domain where gravity darkening has been considered, 
mainly for reproducing the light curves of eclipsing binaries. We may wonder 
if the w-model can be generalized to predict the gravity darkening of a star 
belonging to a binary system. It does but without any (known) analytic solution. 

Let us follow the work of Espinosa Lara & Rieutord (2012). In the radiative 
envelope of an early-type star member of a binary system we can still write the 
conservation of the flux and assume the anti-parallelism of flux and effective 
gravity: 


DivF = 0 and F = - fg oS 
but now the effective gravity comes from the 3D potential: 


GMi 


GM-2 


1 /->2 2 / • 2 /i • 2 

— -ilr (sin Asm ip 


\J a 2 + r 2 — 2ar cos 6 
cos 2 6) - 


Mo 


-frr cos 8 , 


(37) 


Mi + M 2 

where M\ and M 2 are the masses of the two stars, ’a’ is the distance between the 
two stellar centres and 12 is the orbital angular velocity. The orbit is assumed 
circular. Let us write Di v(fg e g) = 0 as 


V • g e ff 

^eff 


n • V In / 


(38) 
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where we set g c g = —g e gn. 

We consider the three-dimensional curve C(9o,<po) that starts at the centre 
of the star with the initial direction given by (0o> and that is tangent to n 
at every point. C(9 0 , ip 0 ) is therefore a field line of the effective gravity field. 
The value of / at a point r along the curve can be calculated as a line integral 


AH = /o' 


C(9 0 ,ip 0 ) 


Despite much efforts no analytical expression could be found for /. Expression 
(39) is thus integrated numerically. 

One interesting result of this approach, is that there is not a one-to-one 
relation between effective gravity and effective temperature. Indeed, because 
of the absence of symmetry of the star (except of the equatorial one if the 
obliquity is zero), two different points of the stellar surface may have the same 
effective gravity but a different effective temperature. This property comes from 


expression (391: the path integrals that lead to two points of identical effective 


gravity are not necessarily the same and can lead to different values of the 
flux. This property is illustrated in Fig. [8j In this figure we see that the curve 
T e « = f(gefi) is not smooth because similar values of g e g lead to different values 
of T e g. Fortunately, these variations are small. 


As in Espinosa Lara & Rieutord (2012), we define q as the mass ratio, and 


evaluate the filling of the Roche lobe by the radius p of the star along the line 
joining the stellar centres, taking the distance between the star centre and the 
Lagrange Li point as unity. Hence, a star filling its Roche lobe has p = 1 while 
the one filling it at 95% has p = 0.95. The different positions where the same 
effective temperature are found, is illustrated in the two cases shown in Fig. [9] 
There we see that the curves of isoflux are not simple curves over the stellar 
surface. 

As shown by Djurasevic et al. |20C)3j, the determination of the /3-exponent 
from the light curves of semi-detached binaries is almost impossible since mag¬ 


netic spots induce similar variations (see Fig. 10 and 11). 


5 Conclusions 

To conclude these notes, I would like to stress a few points on gravity darkening: 

• As far as non-magnetic early-type stars are concerned, gravity darkening has 
no longer to be proved. The use of the /3-model, which is not physically sound 
can be left aside and replaced by the w-model, which has the advantage of 
giving a direct estimate of the w parameter. 

• As far as late-type stars or giant stars are concerned, the situation is much 
more uncertain. The problem is indeed more difficult both on the theoretical 
and observational sides. On the theoretical side, the absence of any univer¬ 
sally accepted model of turbulent rotating convection impedes any serious 
prediction on the latitude dependence of the convective flux. Observational 
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Fig. 8. From Espinosa Lara & Rieutord (2012), correlation of effective temperature 


and effective gravity in the primary early-type star of a binary system. Mass ratio 
is unity and the star fills the Roche lobe at 95% (see text for our definition). The 
correlation may be represented by a /3-exponent of 0.229. The solid line shows a linear 
fit, the dashed line the von Zeipel law, and pluses are from our generalized cj-model. 


•H 



z 



1,68 

1.60 

1.52 

1.36 


1.28 

1.20 

1.12 


Espinosa Lara & Rieutord 


(2012): distribution of the effective tempera- 


Fig. 9. From 

ture at the surface of a tidally distorted star. Left: q = 1 and p = 0.8. Right: q = 1 and 
p = 0.95. 


constraints are therefore most welcome. However, this is not a simple mat¬ 
ter either. Convective envelopes are usually harboring magnetic fields which 
can disturb the flux distribution. Ideally, the surface of these stars should be 
constrained by both interferometers and Zeeman-Doppler Imaging so as to 
disentangle the effects. 

• Finally, for both type of (single) stars, we may recommend the following 
scheme of hypothesis and measurements. First assume the axi- and equatorial 
symmetry of the star. Then, if the star is centrally condensed (like an early- 
type or a giant one), adopt the Roche model. If the star is not centrally 
condensed (like a late-type star of the main sequence), a bipolytropic model 
is fine. Such a model, which fits the radiative and convective zones with a 
polytrope, just depends on three parameters, mass, equatorial radius and 
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TV Cas 




PHASE - 0.25 

Fig. 10. A model of TV Cas that leads to /? = 0.15 from fitting the light curve (from 
Djurasevic et al.|2003 |. 

TV Cas 



PHASE-0.25 


Fig. 11. The second model with a spot and /3 = 0.25 for TV Cas (from Djurasevic 
|et al.|2 003); the difference between the calculated and observed light curve is the same 
as with the model of Fig. mi 


u>, just like the Roche model. Then, the flux or the effective temperature 
T e g(6) can be derived after an expansion on the spherical harmonic basis 
along with an atmosphere model used for the determination of the limb 
darkening effect. The gravity darkening law can then be evaluated from the 
curve (or correlation) T e g(6) versus g e g(6). 
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Appendix: Angular velocity with respect to critical 
rotation 


In this appendix we discuss the correspondance between two definitions of the 
scaled angular velocity. The first one is the one we used in the text, namely 



where is the orbital angular velocity for an orbit at the actual equatorial 
radius of the star. 

The other definition is based on the Roche model and considers the angular 
velocity ft c such that the rotation on the equatorial radius is keplerian. This 
latter definition is a true critical angular velocity, while the previous one is a 
keplerian velocity at the actual equatorial radius. However, the critical angular 
velocity is model dependent, this is why we have to mention Roche’s model 
for the definition of fi c . The first definition does not need any model, but it is 
not the exact critical angular velocity. This latter quantity cannot in general 
be computed a priori with just a given spherical model of a star. It needs a full 
computation of the structure at the actual critical velocity and is thus an output 
of 2D models like ESTER ones (e.g. Espinosa Lara fe Rieutord||2013 ). 

So we now only consider Roche models where all quantities can be derived in a 
simple manner. We recall that the polar R p and equatorial R e of an equipotential 
of a star rotating at angular velocity fl are related by 


= ^ 1 

Rp Re 2 

Then, the critical angular velocity f2 c and the critical equatorial radius R ec 
related by 


(40) 

are 


and hence 


ni 


GM 

R3 


(41) 


Rec = \Rp (42) 

at critical rotation. 

If the rotation is subcritical, the Roche model gives the following relation 
between R e and R p 

Rp = Re (l + y) (43) 

But we may take 12 C as the scale of the rotation rate and set 


w = 


n 

G c 


(44) 
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From the preceding definitions we get the relation between Co and ui, namely 


00 = LO 



-3/2 


(45) 


We note that if w = 1 then Co = 1 as expected. We also note that if 42 is 
subcritical, then R e < 3R p /2 and therefore 42^ > 42 c , which implies that we 
always have 


Co > uj 


(46) 


Equation (45) shows that it is easy to compute Co from oo but the opposite is 
a little more complicated since a cubic equation must be solved. Setting \ = 
arcsinw, we find 


to = 



sin(x/3) - 2 


(47) 
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